T cell

Load in libraries

library(tidyverse)
library(readxl)
library(mgcv)
library(emmeans)
library(ggeffects)
library(lme4)
library(gamlss)
library(hrbrthemes)
library(DT)
library(ggpubr)
library(cowplot)
library(Cairo)
source("scripts/utils/utils.R")
modernacolors()

## <environment: R_GlobalEnv>

Load data

# color panel
panel0 <-
  c(    modernadarkgray,
    modernadarkblue2,
    modernaorange,
    modernapurple,
    modernagreen,
    modernablue,
    modernared
  )
dosepanel <- c("#cb6939", "#8264cb")
d <- read_excel("processed_data/Tcell.xlsx")
dt <- d %>%
  mutate(
    animal_id = factor(animal_id),
    interval = factor(interval),
    day = factor(day, levels = c(64, 87, 143, 227)),
    interval = factor(interval, levels = c("PBS", 0:8))
  )

Model fitting for CD4 data

# data analysis for CD4 cell Th1 responses
dt4 <- dt %>%
  filter(cell_type == "CD4" & interval != "PBS") %>%
  droplevels()

# fit zero inflated beta regression
con1 <- gamlss::gamlss.control(n.cyc = 100, trace = FALSE)
combos <- expand.grid(unique(dt4$cytokine), unique(dt4$stimulation))

# split data by cytokine and stimulation 
datas <-
  map2(combos$Var1, combos$Var2, ~
         dt4 %>% filter(cytokine == .x & stimulation == .y)) %>%
  setNames(paste(combos$Var1, combos$Var2, sep= "_" )) 

# function to fit zero inflated beta regression
fit_fun <- function(dt, compare = interval) {
  data = dt
  if (compare == "interval") {
    formula = as.formula("~ interval | day")
  } else{
    formula = as.formula("~ day|interval")
  }
  tryCatch({
    fit <- gamlss::gamlss(
      prop ~  day * interval,
      sigma.fo = ~ day * interval,
      family = BEZI,
      data = data,
      trace = F,
      control = con1
    )
    em_fit <-
      emmeans::emmeans(ref_grid(fit, cov.keep = c("interval", "day")),
                       specs = formula ,
                       type = "response")
  },
  error = function(e) {
    fit <- gamlss::gamlss(
      prop ~  interval + day,
      sigma.fo = ~ day + interval,
      family = BEZI,
      data = data,
      trace = F,
      control = con1
    )
    em_fit <-
      emmeans(ref_grid(fit, cov.keep = c("interval", "day")),
              specs = formula ,
              type = "response")
  })
  return(list(fit = fit, emm_fit = em_fit))
}

# obtain the model 
set.seed(100)
beta_fits <-
  purrr::map(datas, ~ fit_fun(dt = .x, compare = "interval"))
fits <- purrr::transpose(beta_fits)[[1]] 

emmt_dat2 <- 
  purrr::transpose(beta_fits)[[2]] %>%
  imap(~.x %>% 
         as_tibble() %>%
         mutate(info = .y)) %>%
  bind_rows() %>%
  separate(info, c("cytokine", "coating"), "_")

Residual diagnostics for CD4+ data

assumption_plots <-
  imap(
    fits ,
    ~ assumption_check(.x, .x$y) %>%
      plot_grid(
        ggdraw() + draw_label(
          paste0(str_replace(.y, "_", ", " ), "-specific"),
          x = 0.03,
          hjust = 0
        ),
        .,
        ncol = 1,
        rel_heights = c(0.1, 1),
        align = "v"
      )
  )
ggarrange(plotlist = assumption_plots)

Supplementary Figure 5AB: CD4 fold change and heatmap

#-------------------------------------------------------------------------------
# Supplementary Figure 5A: CD4 fold change
#-------------------------------------------------------------------------------
plots4 <-
  emmt_dat2 %>%
  filter(coating %in% c("S1", "S2")) %>%
   mutate(
    day = forcats::fct_recode(
      factor(day),
      "1 week post-boost" = "64",
      "4 weeks post-boost" = "87",
      "12 weeks post-boost" = "143",
      "24 weeks post-boost" = "227"
    )
  ) %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
     "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  split(.$cytokine) %>%
  imap(
    ~ .x %>%
      filter(day == "24 weeks post-boost") %>%
      ggplot(aes(
        x = interval, y = response * 100, color = interval
      )) +
      geom_jitter(position = position_dodge(width = 0.5), size = 2)  +
      geom_errorbar(
        aes(ymin = asymp.LCL * 100, ymax = asymp.UCL * 100),
        width = .1,
        position = position_dodge(width = 0.5),
        size = 1
      ) +
      scale_color_manual(values = panel0) +
      theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
      facet_grid(day ~ coating, space = "free_x", scales = "free") +
      labs(title = paste0("CD4+ T cells (", .y, ")"),
           x = "", y = "Estimated CD4+ T cells, %") +
      theme(
        axis.text.x = element_text(
          angle = 45,
          hjust = 1,
          vjust = 1
        ),
        legend.position = "None",
        plot.title = element_text( face = "bold"),
        panel.spacing = unit(0, "lines")
      )
  )



set.seed(100)
emmt_dat <-
  purrr::transpose(beta_fits)[[2]] %>%
  imap(
    ~ contrast(
      regrid(regrid(.x, "response"), transform = "log10"),
      "pairwise",
      infer = TRUE,
      adjust = "mvt",
      ratio = FALSE,
      level = .95
    ) %>% as_tibble() %>%
      mutate(info = .y)
  ) %>%
  bind_rows() %>%
  separate(info, c("cytokine", "coating"), "_")

emmt_dat_reverse <-
  emmt_dat %>%
  mutate("significant or not" = ifelse(asymp.LCL < 0 &
                                         asymp.UCL > 0, "no", "yes")) %>%
  dplyr::select(contrast, day, estimate, cytokine, coating, `significant or not`) %>%
  separate(contrast, c("g2", "g1"), " - ") %>%
  mutate(estimate = -estimate)

# heat map-----------
heatmaps4 <-
  emmt_dat %>%
  mutate("significant or not" = ifelse(asymp.LCL < 0 &
                                         asymp.UCL > 0, "no", "yes")) %>%
  dplyr::select(contrast, day, estimate, cytokine, coating, `significant or not`) %>%
  separate(contrast, c("g1", "g2"), " - ") %>%
  bind_rows(emmt_dat_reverse) %>%
  mutate(g1 = str_remove(g1, "interval"),
         g2 = str_remove(g2, "interval")) %>%
  filter(coating %in% c("S1", "S2")) %>%
  split(.$cytokine) %>%
  imap(
    ~ complete(.x, g2, g1, day, cytokine, coating) %>%
      mutate("significant or not" = ifelse(
        is.na(`significant or not`), "no", `significant or not`
      )) %>%
      dplyr::rename(`Estimated fold change` = `estimate`) %>%
      filter(day == 227) %>%
      ggplot(aes(
        x = g2, y = g1, fill = `Estimated fold change`
      )) +
      scale_x_discrete(position = "top") +
      scale_fill_gradient2(
        na.value = "grey25",
        breaks = log10(c(1 / 10, 1 / 3,  1, 3, 10)),
        labels =  c("1/10", "1/3", 1, 3, 10),
        high = scales::muted("red"),
        low = scales::muted("blue")
      ) +
      theme_ipsum(base_family="sans") +
      geom_tile() +
      geom_point(
        aes(shape = `significant or not`),
        color = "black",
        size = 2
      ) +
      scale_shape_manual(values = c(NA, 8)) +
      labs(shape = "significantly\ndifferent?", x = "Comparison", y = "",
            #title = paste("CD4+ T cells (", .y, ")")
           ) +
      facet_wrap(
        ~ coating,
        scales = "free",
        strip.position = "top",
        nrow = 1
      ) +
      theme(
        strip.placement = "outside"
      )
  )

Supplementary Table 6. Statistical comparisons of percentage of SARS-CoV-2 spike-specific IFNγ, IL-2, and TNFα producing CD4+ T cells through 24 weeks following dose 2.

# CD 4 table 6
set.seed(100)
emmt_table <-
  purrr::transpose(beta_fits)[[2]] %>%
  imap(
    ~ contrast(
      regrid(regrid(.x, "response"), transform = "log10"),
      "pairwise",
      infer = TRUE,
      adjust = "mvt",
      level = .95
    ) %>% as_tibble() %>%
      mutate(info = .y)
  ) %>%
  bind_rows() %>%
  separate(info, c("cytokine", "coating"), "_")

Tcell_CD4_ST6 <-
  emmt_table %>%
  as_tibble() %>%
  dplyr::select(contrast, ratio, day, cytokine, coating, p.value) %>%
  mutate(contrast = str_remove_all(contrast, "interval")) %>%
  mutate(
    day = forcats::fct_recode(
      factor(day),
      "7 (1)" = "64",
      "30 (4)" = "87",
      "86 (12)" = "143",
      "169 (24)" = "227"
    )
  ) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "ratio",
    "Cytokine" = "cytokine",
    "Coating" = "coating",
    "Adjusted P-value" = "p.value"
  ) %>%
  mutate(
    `Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
    `Fold change` = round(`Fold change`, digits = 6)
  ) %>%
  filter(`Adjusted P-value` < 0.05)

Tcell_CD4_ST6 %>%
  DT::datatable(
    caption = "Statistical comparisons of percentage of SARS-CoV-2
                spike-specific IFNγ, IL-2, and TNFα producing CD4+ T cells through 24 weeks
                following dose 2"
  )

Supplementary Figure 12 Statistical comparisons of percentage of SARS-CoV-2 spike-specific IFNγ, IL-2, and TNFα producing CD4+ T cells through 24 weeks following dose 2 (Supplementary Table 6 visualization)

Each dot corresponds to the estimated fold change for the contrasted prime-boost dosing interval. Error bars are multivariate t adjusted 95% CIs on the estimated contrast.

cd4 <-
  emmt_table %>%
  as_tibble() %>%
  separate(contrast, c("a", "b"), sep = " / ") %>%
  mutate(a = str_remove_all(a , "interval"),
         b = str_remove_all(b , "interval")) %>%
  mutate(
    a = case_when(
      str_detect(a, "1") ~ paste(a, "week"),
      str_detect(a, "0") ~ str_replace(a, "0", "Prime only"),
      TRUE ~ paste(a, "weeks")
    ),
    b = case_when(
      str_detect(b, "1") ~ paste(b, "week"),
      str_detect(b, "0") ~ str_replace(b, "0", "Prime only"),
      TRUE ~ paste(b, "weeks")
    )
  ) %>%
  mutate(contrast = paste(a, "/", b)) %>%
  mutate(day = case_when(
    day == "64" ~ "7_1",
    day == "87" ~ "30_4",
    day == "143" ~ "86_12",
    day == "227" ~ "169_24"
  )) %>%
  mutate(contrast = str_remove_all(contrast, "interval")) %>%
  filter(coating == "S1" | coating == "S2") %>%
  mutate(group = paste(coating, "specific", cytokine)) %>%
  mutate(
    "Statistically Significant" = ifelse(p.value < 0.05, "yes", "no"),
    `Statistically Significant` = factor(`Statistically Significant`, levels = c("yes", "no"))
  ) %>%
  filter(!is.na(`Statistically Significant`)) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "ratio",
    "Adjusted P-value" = "p.value"
  ) %>%
  mutate(group = factor(
    group,
    levels = c(
      "S1 specific IFNg",
      "S1 specific IL-2",
      "S1 specific TNFa",
      "S2 specific IFNg",
      "S2 specific IL-2",
      "S2 specific TNFa"
    ),
    ordered = TRUE,
    labels = c(
      expression(S1 ~ specific ~ IFN ~ gamma),
      expression(S1 ~ specific ~ IL - 2),
      expression(S1 ~ specific ~ TNF ~ alpha),
      expression(S2 ~ specific ~ IFN ~ gamma),
      expression(S2 ~ specific ~ IL - 2),
      expression(S2 ~ specific ~ TNF ~ alpha)
    )
  ))


cd4_emmeans <-
  cd4 %>%
  split(.$day) %>%
  imap(
    ~ .x %>%
      ggplot(
        aes(x = `Pairwise comparison` , y = `Fold change`, color = `Statistically Significant`)
      ) +
      geom_point(size = 3)  +
      geom_errorbar(
        aes(ymin = asymp.LCL, ymax = asymp.UCL),
        width = .1,
        size = 1
      ) +
      theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
      scale_y_continuous(
        breaks = (c(1, 3, 10)),
        labels = c(1, 3, 10),
        trans = "sqrt"
      ) +
      geom_hline(
        yintercept = 1,
        color = "darkcyan",
        lty = "dashed"
      ) +
      labs(
        y = "Estimated Fold Change",
        x = "Prime-boost dosing interval",
        title = paste0("Approximate ",
                       str_split(.y, "_")[[1]][2],
                       " week post-dose 2")
      ) +
      facet_grid(
        ~ group,
        scales = "free",
        space = "free",
        labeller = label_parsed
      ) +
      coord_flip() +
      theme(
        strip.text = element_text(size = 12),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        legend.position = "top"
      )
    
  )

legend4 <- get_legend(cd4_emmeans$`7_1`)

pgrid4 <-
  cowplot::plot_grid(
    cd4_emmeans$`7_1` + theme(legend.position = "none"),
    cd4_emmeans$`30_4` + theme(legend.position = "none"),
    cd4_emmeans$`86_12` + theme(legend.position = "none"),
    cd4_emmeans$`169_24` + theme(legend.position = "none"),
    ncol = 1,
    labels = c('A', 'B', 'C', 'D'),
    align = 'hv',
    rel_heights = c(1, 1, 1, 1)
  )

p4 <- plot_grid(legend4,
                pgrid4,
                ncol = 1,
                rel_heights  = c(.05, 1))
p4

ggsave(
  filename = "figures_tables/Tcell/cd4_emmean.pdf",
  plot = p4,
  width = 26,
  height = 24,
  device = cairo_pdf
)

ggsave(
  filename = "figures_tables/Tcell/cd4_emmeans.png",
  plot = p4,
  width = 26,
  height = 24
)

Model fitting for CD8 data

# data analysis for CD8 cell Th1 responses
dt8 <- subset(dt, cell_type == "CD8" & interval !="PBS") %>%
  droplevels()

# fit zero inflated beta regression
con1 <- gamlss.control(n.cyc=40, trace = FALSE)
combos <- expand.grid(unique(dt8$cytokine), unique(dt8$stimulation))

# split data by cytokine and stimulation 
datas8 <-
  map2(combos$Var1, combos$Var2, ~
         dt8 %>% filter(cytokine == .x & stimulation == .y)) %>%
  setNames(paste(combos$Var1, combos$Var2, sep= "_" )) 

# obtain the model 
set.seed(100)
beta_fits8 <-
  purrr::map(datas8, ~ fit_fun(dt = .x, compare = "interval"))

# obtain emmeans
emmt_ind8 <- 
  purrr::transpose(beta_fits8)[[2]]%>%
  imap(~.x %>% 
         as_tibble()%>%
         mutate(info = .y)) %>%
  bind_rows()%>%
  separate(info, c("cytokine", "coating"), "_")

Residual diagnostics for CD8+ data

fits <- purrr::transpose(beta_fits8)[[1]] 

assumption_plots <-
  imap(
    fits ,
    ~ assumption_check(.x, .x$y) %>%
      plot_grid(
        ggdraw() + draw_label(
          paste0(str_replace(.y, "_", ", " ), "-specific"),
          x = 0.03,
          hjust = 0
        ),
        .,
        ncol = 1,
        rel_heights = c(0.1, 1),
        align = "v"
      )
  )
ggarrange(plotlist = assumption_plots)

Supplementary Figure 5CD: CD8 fold change

set.seed(100)
emmt_dat8 <- 
  purrr::transpose(beta_fits8)[[2]]%>%
  imap(~contrast(regrid(regrid(.x, "response"), transform = "log10"), 
                 "pairwise",
                 infer = TRUE,
                 adjust = "mvt",
                 ratio = FALSE,
                 level = .95)%>% as_tibble() %>%
         mutate(info = .y)) %>%
  bind_rows()%>%
  separate(info, c("cytokine", "coating"), "_")

emmt_dat_reverse8 <- 
  emmt_dat8 %>%
  mutate("significant or not" = ifelse(asymp.LCL<0 & asymp.UCL >0, "no", "yes"))%>%
  dplyr::select(contrast, day, estimate, cytokine, coating, `significant or not`)%>%
  separate(contrast, c("g2", "g1"), " - ")%>%
  mutate(estimate = -estimate)

plots8 <-
  emmt_ind8 %>%
   mutate(
    day = forcats::fct_recode(
      factor(day),
      "1 week post-boost" = "64",
      "4 weeks post-boost" = "87",
      "12 weeks post-boost" = "143",
      "24 weeks post-boost" = "227"
    )
  ) %>%
  mutate(
    interval = forcats::fct_recode(
      interval,
      "8-week interval" = "8",
      "6-week interval" = "6",
      "4-week interval" = "4",
      "3-week interval" = "3",
      "2-week interval" = "2",
      "1-week interval" = "1",
      "Prime only" = "0"
    )
  ) %>%
  filter(coating %in% c("S1", "S2")) %>%
  split(.$cytokine) %>%
  imap(
    ~ .x %>%
      filter(day == "24 weeks post-boost") %>%
      ggplot(aes(
        x = interval, y = response * 100, color = interval
      )) +
      geom_jitter(position = position_dodge(width = 0.5), size = 2)  +
      geom_errorbar(
        aes(ymin = asymp.LCL * 100, ymax = asymp.UCL * 100),
        width = .1,
        position = position_dodge(width = 0.5),
        size = 1
      ) +
      scale_color_manual(values = panel0) +
      theme_bw() +  theme(panel.grid.major = element_line(colour = "grey100")) +
      facet_grid(day ~ coating, space = "free_x", scales = "free") +
      labs(
        title = paste("CD8+ T cells (", .y, ")"),
        x = "",
        y = "Estimated CD8+ T cells, %"
      ) +
      theme(
        axis.text = element_text(size = 9),
        axis.text.x = element_text(
          angle = 45,
          hjust = 1,
          vjust = 1
        ),
        plot.title = element_text( face = "bold"),
        panel.spacing = unit(0, "lines"),
        legend.position = "None"
      )
    
  )


### Supplementary Figure 5D: CD8 heatmap
heatmaps8 <-
  emmt_dat8 %>%
  mutate("significant or not" = ifelse(asymp.LCL < 0 &
                                         asymp.UCL > 0, "no", "yes")) %>%
  dplyr::select(contrast, day, estimate, cytokine, coating, `significant or not`) %>%
  separate(contrast, c("g1", "g2"), " - ") %>%
  bind_rows(emmt_dat_reverse8) %>%
  mutate(g1 = str_remove(g1, "interval"),
         g2 = str_remove(g2, "interval")) %>%
  filter(coating %in% c("S1", "S2")) %>%
  split(.$cytokine) %>%
  imap(
    ~ complete(.x, g2, g1, day, cytokine, coating) %>%
      mutate("significant or not" = ifelse(
        is.na(`significant or not`), "no", `significant or not`
      )) %>%
      dplyr::rename(`Estimated fold change` = `estimate`) %>%
      filter(day == 227) %>%
      ggplot(aes(
        x = g2, y = g1, fill = `Estimated fold change`
      )) +
      scale_x_discrete(position = "top") +
      scale_fill_gradient2(
        na.value = "grey25",
        breaks = log10(c(1 / 10, 1 / 3, 1 / 2, 1, 2, 3, 10)),
        labels =  c("1/10", "1/3", "1/2", 1, 2, 3, 10),
        high = scales::muted("red"),
        low = scales::muted("blue")
      ) +
      theme_ipsum(base_family="sans") +
      geom_tile() +
      geom_point(
        aes(shape = `significant or not`),
        color = "black",
        size = 2
      ) +
      scale_shape_manual(values = c(NA, 8)) +
      labs(shape = "significantly\ndifferent?", x = "Comparison", y = "",
            #title = paste("CD8+ T cells (", .y, ")")
           ) +
      facet_wrap(
        ~ coating,
        scales = "free",
        strip.position = "top",
        nrow = 1
      ) +
      theme(
        strip.placement = "outside"
      )
  )

Each dot in the left panel of the following figure corresponds to the estimated mean for the prime-boost dosing interval. Error bars are 95% CIs on the estimated mean. Right panel shows statistical comparisons between the dosing interval groups, with significant differences denoted by an asterisk if the P-value was less than 0.05.

hmaps <-
  cowplot::plot_grid(
    plotlist = list(heatmaps4[[1]], NULL, heatmaps8[[1]], NULL),
    rel_heights = c(50, 1, 50, 7),
    nrow = 4,
    ncol = 1
  )
lineplots <-
  cowplot::plot_grid(
    plotlist = list(NULL, plots4[[1]],  plots8[[1]]),
    nrow = 3,
    ncol = 1,
    rel_heights = c(2, 10, 10)
  )
figure5 <-
  cowplot::plot_grid(
    plotlist = list(lineplots, hmaps),
    nrow = 1,
    ncol = 2,
    rel_widths = c(0.9, 1.4)
  )
figure5

Supplementary Table 7. Statistical comparisons of percentage of SARS-CoV-2 spike-specific IFNγ, IL-2, and TNFα producing CD8+ T cells through 24 weeks following dose 2

## CD 8 table 7
set.seed(100)
emmt_table8 <-
  purrr::transpose(beta_fits8)[[2]] %>%
  imap(
    ~ contrast(
      regrid(regrid(.x, "response"), transform = "log10"),
      "pairwise",
      infer = TRUE,
      adjust = "mvt",
      level = .95
    ) %>% as_tibble() %>%
      mutate(info = .y)
  ) %>%
  bind_rows() %>%
  separate(info, c("cytokine", "coating"), "_")

Tcell_CD8_ST7 <-
  emmt_table8 %>%
  as_tibble() %>%
  dplyr::select(contrast, ratio, day, cytokine, coating, p.value) %>%
  mutate(contrast = str_remove_all(contrast, "interval")) %>%
  mutate(
    day = forcats::fct_recode(
      factor(day),
      "7 (1)" = "64",
      "30 (4)" = "87",
      "86 (12)" = "143",
      "169 (24)" = "227"
    )
  ) %>%
  filter(p.value < 0.05) %>%
  rename(
    "Pairwise comparison" = contrast,
    "Fold change" = "ratio",
    "Cytokine" = "cytokine",
    "Coating" = "coating",
    "Adjusted P-value" = "p.value"
  ) %>%
  mutate(
    `Adjusted P-value` = round(`Adjusted P-value`, digits = 6),
    `Fold change` = round(`Fold change`, digits = 6)
  ) %>%
  filter(`Adjusted P-value` < 0.05)

Tcell_CD8_ST7 %>%
  DT::datatable(caption = "Statistical comparisons of percentage of SARS-CoV-2 spike-specific IFNγ, IL-2, and TNFα producing CD8+ T cells through 24 weeks following dose 2
")

Supplementary Figure 13 Statistical comparisons of percentage of SARS-CoV-2 spike-specific IFNγ, IL-2, and TNFα producing CD8+ T cells through 24 weeks following dose 2 (Supplementary Table 7 visualization)